#include "CFraction.h"

CFraction::CFraction(shared_ptr<TaylorSeries> pSeries, double c0) : Approximant(pSeries), _c0(c0), _pprogresiveQD(_pSeries, _c0), _degree(0), _x(1.0)
{
}

CFraction::~CFraction(void)
{
}

std::pair<double, double> CFraction::operator()()
{
	auto coeffs = std::make_pair(a(_counter), b(_counter));
	_counter++;
	return coeffs;
}

double CFraction::a( unsigned int n )
{
	return (*_pnumeratorCoeffs)[n] * _x;
}

double CFraction::b( unsigned int n )
{
	if(n == 0)
		return 0.0;
	return 1;
}

double CFraction::evaluate( double x )
{
	return forwardRecurrence(x); // fast but may become unstable
	//return backwardRecurrence(x); // stable but slow
	//return lentz(x); // stable but even slower
}

void CFraction::build( unsigned int n )
{
	_degree = n;
	
	_pnumeratorCoeffs = _pprogresiveQD.partialNumerators(n);	
}

double CFraction::lentz( double x )
{
	_x = x;
	_counter = 0;
	boost::uintmax_t maxTerm = (_degree - 1);
	return boost::math::tools::continued_fraction_a<CFraction>(*this, 0.0, maxTerm);
}

double CFraction::backwardRecurrence( double x )
{
	double c = 1.0;

	std::vector<double>& a = *_pnumeratorCoeffs;

	for(int k = _degree - 1; k > 0; k--)
	{
		c = 1.0 + ((a[k]*x)/c);
	}

	return ( a[0]*x ) / c;
}

double CFraction::forwardRecurrence( double x )
{
	return _numerator.evaluate(x)/_denominator.evaluate(x);
}

// Since on some (actually most) architectures multiply is much faster than divide
// it is useful to convert the convergent into a rational function. This is achieved
// through the forward recurrence method, note however the algorithm may be unstable
// in some cases. 
void CFraction::convert()
{
	int n = degree();
	std::vector<Polynomial> A(n+2), B(n+2);
	
	//initialize
	_x = 1.0;
	A[0].setCoefficient(0,1.0);
	A[1].setCoefficient(0,b(0));
	B[0].setCoefficient(0,0.0);
	B[1].setCoefficient(0,1.0);

	for(int j = 1; j <= n; j++)
	{
		Polynomial aj;
		aj.setCoefficient(1,a(j-1));
		A[j+1] = (A[j]*b(j))+(A[j-1]*aj);
		B[j+1] = (B[j]*b(j))+(B[j-1]*aj);
	}
	_numerator = A[n+1];
	_denominator = B[n+1];

}

